home *** CD-ROM | disk | FTP | other *** search
/ APDL Eductation Resources / APDL Eductation Resources.iso / programs / astronomy / skyview / !SkyView / c / RADec < prev    next >
Encoding:
Text File  |  1993-08-26  |  9.4 KB  |  294 lines
  1. /********************************************************/
  2. /*                         RADec                        */
  3. /* Utility functions for objects with known RA and Dec. */
  4. /********************************************************/
  5.  
  6. #include <math.h>
  7. #include <stdio.h>
  8.  
  9. #include "os.h"
  10. #include "coords.h"
  11. #include "menu.h"
  12. #include "sv_header.h"
  13. #include "radec.h"
  14. #include "datime.h"
  15.  
  16. /* Set declination limit, above which object is consid- */
  17. /* ered so close to the celestial pole that it doesn't  */
  18. /* move during the day:                                 */
  19. #define OBJECT_NOMOVE (REAL)1.57
  20.  
  21. /*------------------ Private functions -----------------*/
  22. static void phenom_details(
  23.             REAL ra, REAL dec, observerstr *ob_ptr, REAL hourang,
  24.             double *cosazptr, int *hourptr, int *minptr);
  25.  
  26. /*------------------- Global variables -----------------*/
  27. static double horalt;  /* 'Altitude' of horizon.        */
  28.  
  29.  
  30. /*--------------------- Functions ----------------------*/
  31.  
  32. /*------------------------------------------------------*/
  33. /*  Function to set altitude which object must exceed   */
  34. /*          in order to be regarded as visible.         */
  35. /*------------------------------------------------------*/
  36.   void radec_sethoriz(REAL altitude)
  37. {
  38.   horalt = (double)altitude;
  39.   return;
  40. }
  41.  
  42. /*------------------------------------------------------*/
  43. /*  Convert from RA & Dec to Alt & Azim following the   */
  44. /*  method of Duffett-Smith p39.  Angles in radians.    */
  45. /*------------------------------------------------------*/
  46.     void radec_altaz(REAL ra, REAL dec, observerstr *ob_ptr, REAL sid,
  47.                      REAL *altptr, REAL *azimptr)
  48. {
  49.   static double altmax = 1.5700;
  50.   static double h;
  51.   static double sind, cosd, sinlat, coslat;
  52.   static double sinalt, alt;
  53.   static double cosazim, azim;
  54.  
  55. /* Calculate local hour angle of object.                */
  56.   h = (double)sid - (double)ra;
  57.  
  58. /* Calculate sin and cos of object's declination and    */
  59. /* observer's latitude.                                 */
  60.   sind   = sin((double)dec);
  61.   cosd   = cos((double)dec);
  62.   sinlat = sin((double)ob_ptr->latit);
  63.   coslat = cos((double)ob_ptr->latit);
  64.  
  65. /* Calculate altitude of object.                        */
  66.   sinalt = sind*sinlat + cosd*coslat*cos(h);
  67.   if (sinalt < -1.0) sinalt = -1.0;
  68.   if (sinalt >  1.0) sinalt =  1.0;
  69.   alt = asin(sinalt);
  70.  
  71. /* Calculate azimuth of object.                         */
  72. /* Special case when object is very nearly overhead.    */
  73.   if (alt > altmax)
  74.   /* Arbitrarily set azimuth to zero.                   */
  75.     azim = 0.0;
  76.   else
  77.   {
  78.     cosazim = (sind - sinlat*sinalt)/(coslat*cos(alt));
  79.     if (cosazim < -1.0) cosazim = -1.0;
  80.     if (cosazim >  1.0) cosazim =  1.0;
  81.     azim = acos(cosazim);
  82.   }
  83.  
  84.   *altptr = (REAL)alt;
  85.   *azimptr= (sin(h) < 0.0 ? (REAL)azim: (REAL)2.0*PI - (REAL)azim);
  86.  
  87.   return;
  88. }
  89.  
  90.  
  91. /*------------------------------------------------------*/
  92. /*  Function to calculate whether the selected object   */
  93. /*  rises, sets or culminates today. Angles in radians. */
  94. /*------------------------------------------------------*/
  95.     void radec_phenomena(REAL ra, REAL dec, observerstr *ob_ptr,
  96.                          double *chp, BOOL *riset_ptr, BOOL *cul_ptr)
  97. {
  98. /* Max zenith distance for object to be visible at      */
  99. /* culmination:                                         */
  100.   REAL zd_vis = PIby2 - (REAL)horalt;
  101.  
  102.   double d_latit  = (double)ob_ptr->latit;
  103.   double d_dec    = (double)dec;
  104.  
  105.  
  106. /* If object is at or very near the celestial pole, it  */
  107. /* does not move, and therefore does not rise, set or   */
  108. /* culminate.                                           */
  109.   if (dec >= OBJECT_NOMOVE  || \
  110.       dec <=-OBJECT_NOMOVE)
  111.   {
  112.     *riset_ptr = FALSE;
  113.     *cul_ptr   = FALSE;
  114.     return;
  115.   }
  116.  
  117. /*              *** Rising and Setting ***              */
  118.  
  119. /*    Object is deemed to 'rise' & 'set' when at an     */
  120. /*         altitude of 'horalt' radians.                */
  121.  
  122. /* Calculate cos of hour angle as object crosses this   */
  123. /* altitude.                                            */
  124.   *chp = (sin(horalt) - sin(d_latit)*sin(d_dec)) / \
  125.                          (cos(d_latit)*cos(d_dec));
  126.  
  127. /* If *chp is not a valid cosine, object does not rise  */
  128. /* or set.                                              */
  129.   if (*chp >=  1.0 || *chp <= -1.0)
  130.     *riset_ptr = FALSE;
  131.   else
  132.     *riset_ptr = TRUE;
  133.  
  134.  
  135. /*                 *** Culmination ***                  */
  136.  
  137. /* Object is visible at culmination if less than zd_vis */
  138. /* radians from the zenith when local hour angle of     */
  139. /* Aries = rt ascension of object.                      */
  140.   *cul_ptr = \
  141.     dec > ob_ptr->latit - zd_vis && \
  142.     dec < ob_ptr->latit + zd_vis;
  143.  
  144.   return;
  145.  
  146. }
  147.  
  148.  
  149. /*------------------------------------------------------*/
  150. /*    Function to calculate details of culmination.     */
  151. /*               All angles in radians.                 */
  152. /*------------------------------------------------------*/
  153.     void radec_cul_details(REAL ra, REAL dec, observerstr *ob_ptr,
  154.                            REAL *al, REAL *az, int *hourptr, int *minptr)
  155. {
  156.   REAL angdiff = dec - ob_ptr->latit;
  157.  
  158. /* Direction to look to see object at culmination.      */
  159.   if (angdiff >= (REAL)0.0)
  160.   {
  161.     *al = PIby2 - angdiff;
  162.     *az = (REAL)0.0;
  163.   }
  164.   else
  165.   {
  166.     *al = PIby2 + angdiff;
  167.     *az = PI;
  168.   }
  169.  
  170. /*                 Time of culmination.                 */
  171. /* Object culminates when local siderial time = ra.     */
  172.  
  173. /* Calculate (one) civil time today at which the local  */
  174. /* siderial time has the required value.                */
  175.   datime_time_today(ra, ob_ptr, hourptr, minptr);
  176.  
  177.   return;
  178. }
  179.  
  180.  
  181. /*------------------------------------------------------*/
  182. /*      Function to calculate details of rising.        */
  183. /*               All angles in radians.                 */
  184. /*------------------------------------------------------*/
  185.     void radec_rise_details(REAL ra, REAL dec, observerstr *ob_ptr,
  186.                             double ch, REAL *az, int *hourptr, int *minptr)
  187. {
  188.   double cosazim;      /* Cosine of azimuth of rising.  */
  189.  
  190. /* Get hour angle of object on rising.                  */
  191.   REAL hourang = -(REAL)acos(ch);
  192.  
  193. /* Get hour, min and azimuth of rising.                 */
  194.   phenom_details(ra, dec, ob_ptr, hourang, &cosazim, hourptr, minptr);
  195.   *az = (REAL)acos(cosazim);
  196.  
  197.   return;
  198. }
  199.  
  200.  
  201. /*------------------------------------------------------*/
  202. /*      Function to calculate details of setting.       */
  203. /*               All angles in radians.                 */
  204. /*------------------------------------------------------*/
  205.     void radec_set_details(REAL ra, REAL dec, observerstr *ob_ptr,
  206.                            double ch, REAL *az, int *hourptr, int *minptr)
  207. {
  208.   double cosazim;  /* Cosine of azimuth of setting.     */
  209.  
  210. /* Get hour angle of object on setting.                 */
  211.   REAL hourang = (REAL)acos(ch);
  212.  
  213. /* Get hour, min and azimuth of setting.                */
  214.   phenom_details(ra, dec, ob_ptr, hourang, &cosazim, hourptr, minptr);
  215.   *az = TWO_PI - (REAL)acos(cosazim);
  216.  
  217.   return;
  218. }
  219.  
  220. /*------------------------------------------------------*/
  221. /* Function to calculate details of rising and setting. */
  222. /* Based on function in Duffett-Smith, p106.            */
  223. /*------------------------------------------------------*/
  224. static void phenom_details(
  225.             REAL ra, REAL dec, observerstr *ob_ptr, REAL hourang,
  226.             double *cosazptr, int *hourptr, int *minptr)
  227. {
  228.   double d_latit  = (double)ob_ptr->latit;
  229.   double d_dec    = (double)dec;
  230.   REAL lst;        /* Local siderial time of phenomenon.*/
  231.  
  232. /* Get local siderial time of phenomenon.               */
  233.   lst = ra + hourang;
  234.  
  235. /* Find civil time corresponding to this siderial time. */
  236.   datime_time_today(lst, ob_ptr, hourptr, minptr);
  237.  
  238. /* Find cosine of azimuth of phenomenon.                */
  239.   *cosazptr = (sin(d_dec) - sin(d_latit)*sin(horalt)) / \
  240.                            (cos(d_latit)*cos(horalt));
  241.  
  242. /* Make sure *cosazptr is a valid cosine.               */
  243.   if (*cosazptr >  1.0) *cosazptr =  1.0;
  244.   if (*cosazptr < -1.0) *cosazptr = -1.0;
  245.  
  246.   return;
  247. }
  248.  
  249.  
  250. /*------------------------------------------------------*/
  251. /*   Function to build a string giving RA and Dec in    */
  252. /*          text form in conventional units.            */
  253. /* String pointed to by textptr is assumed long enough. */
  254. /*------------------------------------------------------*/
  255.   void radec_text(REAL ra, REAL dec, char *textptr)
  256. {
  257.   REAL ram;    /* Right Ascension in minutes.          */
  258.   REAL decm;   /* Declination in minutes.              */
  259.   int  hrai;   /* Hours of RA (integer).               */
  260.   int  mrai;   /* Minutes of RA (integer).             */
  261.   int  ddeci;  /* Degrees of Dec (integer).            */
  262.   int  mdeci;  /* Minutes of Dec (integer).            */
  263.   int ns;      /* North or South declination.          */
  264.  
  265. /* Get RA in conventional units.  Round to nearest min.*/
  266.   ram  = ra * (REAL)720.0 / PI + (REAL)0.5;
  267.   if (ram < ZERO) ram = ZERO;
  268.   hrai = (int)ram / 60;
  269.   mrai = (int)ram % 60;
  270.   if (hrai > 23) {hrai = 0;
  271.                   mrai = 0;}
  272.  
  273. /* Get Dec in conventional units. Round to nearest min.*/
  274.   decm = dec * (REAL)60.0 / CONV;
  275.   if (decm < ZERO)
  276.   {
  277.     decm = -decm;
  278.     if (decm < ZERO) decm = ZERO;
  279.     ns = 'S';
  280.   }
  281.   else
  282.     ns = 'N';
  283.   decm += (REAL)0.5;
  284.   ddeci = (int)decm / 60;
  285.   mdeci = (int)decm % 60;
  286.   if (ddeci > 89) {ddeci = 90;
  287.                    mdeci =  0;}
  288.  
  289.   sprintf(textptr, "RA %02ih %02im   Dec %02i° %02i' %c",
  290.           hrai, mrai, ddeci, mdeci, ns);
  291.  
  292.   return;
  293. }
  294.